library(SoupX)
library(Seurat)
library(tidyverse)

browseVignettes('SoupX')
# quick start in my dataset -------

# we need a 10X outs dir containing raw_ and filtered_ matrices
sc <- load10X()

sc <- autoEstCont(sc)

# if not set roundToInt will produce float number by default
out <- adjustCounts(sc, roundToInt = TRUE)

# in detail -------
tmpDir = tempdir(check = TRUE)

# table of droplets / table of counts
download.file("https://cf.10xgenomics.com/samples/cell-exp/2.1.0/pbmc4k/pbmc4k_raw_gene_bc_matrices.tar.gz", 
              destfile = file.path(tmpDir, "tod.tar.gz"))
download.file("https://cf.10xgenomics.com/samples/cell-exp/2.1.0/pbmc4k/pbmc4k_filtered_gene_bc_matrices.tar.gz", 
              destfile = file.path(tmpDir, "toc.tar.gz"))
untar(file.path(tmpDir, "tod.tar.gz"), exdir = tmpDir)
untar(file.path(tmpDir, "toc.tar.gz"), exdir = tmpDir)

sc <- load10X(tmpDir)

# if you only have filtered matrix (most database)
toc <- sc$toc

# Calculate soup profile
soupProf <- data.frame(
  row.names = rownames(toc),
  est = rowSums(toc)/sum(toc),
  counts = rowSums(toc))

scNoDrops <- toc %>%
  SoupChannel(toc, calcSoupProfile = FALSE) %>%
  setSoupProfile(soupProfile = soupProf)

scNoDrops %>%
  setClusters(setNames(PBMC_metaData$Cluster,
                       rownames(PBMC_metaData))) %>%
  autoEstCont()

# add clustering metadata ---------
data(PBMC_metaData)
sc <- setClusters(sc,
                  setNames(PBMC_metaData$Cluster,
                           rownames(PBMC_metaData)))

sc <- setDR(sc, PBMC_metaData[colnames(sc$toc), c("RD1", "RD2")])

# visual sanity check ---------
ggplot(PBMC_metaData, aes(RD1, RD2)) +
  geom_point(aes(colour = Annotation), size = 0.2) +
  ggtitle("PBMC 4k Annotation") + 
  guides(colour = guide_legend(override.aes = list(size = 1)))

PBMC_metaData$IGKC <- toc['IGKC',]

ggplot(PBMC_metaData, aes(RD1, RD2)) +
  geom_point(aes(colour = IGKC > 0))

# manually set contamination fraction ------
# set a high threshold will remove some true counts
# but will retain good marker genes
# 20% threshold will remove >99% of soup
sc <- setContaminationFraction(sc, 0.2)

# auto-set contamination fraction
sc <- autoEstCont(sc)

# correcting expression profile ------
out <- adjustCounts(sc, roundToInt = TRUE)

plotChangeMap(sc, out, "IGKC")
